Toy model for the mean-field theory of hard-sphere liquids 
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I. INTRODUCTION 



The theory of classical liquids [QJ^] received recently 
an important stimulus from the the theory of structural 
glasses pH-p"l]| . In the pioneering series of papers by Kirk- 
patrick and Thirumalai |lj,[ll| the possible connection of 
structural-glass transition with the spin-glass transition 
in p-spin models was put forward. The analogy was then 
developed e. g. in the problem of minimally correlated 
sequences, which was shown to possess a glassy behavior 
without quenched randomness | l^ , |l5| . 

However, it would be desirable to put this analogy 
further. One of the difficulties occurring in structural 
glasses, when compared to spin glasses, comes from the 
absence of any kind of analytically solvable mean-field 
version of the model. In the case of spin glass, the role is 
played by the fully-connected Ising model, which is solv- 
able very easily. The disordered version is the well-known 
Sherrington-Kikpatrick model |l6[, w hose understanding 
is now very close to be complete ]17|] . 

On the contrary, no analytical solution of a "mean- 
field" liquid is known, as far as we know. It can even 
sound not very reasonable to speak about a mean-field 
liquid, because the relevant high-density phase is char- 
acterized by strong short-range correlations, which can 
hardly be replaced by an effective medium. So, the mean- 
ing of the mean field should be better specified. In our 
investigation, we will understand by mean-field the sit- 
uation which occurs in very high dimension, D — > oo. 
The purpose of the present work is to introduce a sim- 
ple model of a liquid, which is analytically solvable in 
the limit of infinite dimension, at least in a certain well- 
defined range of densities. 



II. SIMPLIFIED HNC EQUATIONS 

We consider a liquid composed of hard spheres with the 
diameter 1. There is only one independent state variable, 
which is the spatial density of particles p. 



The configuration of the liquid is described by the ra- 
dial pair distribution function g(r) — h(r ) + 1. In the 
hypernetted chain (HNC) approximation [|18| we have a 
closed set of equations for the correlation function h(r) 



h(r) + 1 = exp(W(r) - pu{r)) 



(1) 



W(p) = 



ph\p) 
1 + ph{ P ) 



The potential is U(r) = for r > I and U(r) = oo for 
r < 1. These equations can be interpreted as conditions 
for minimization of the free energy functional Q 

T[h] = 

p 2 [ dr r D -\{h(r) + l)(ln(/j(r) +!)-! + U(r)) + 1) (2) 



(27T) 



dpp^L^phip)) 



The function L 3 (x) 
» oo for x — > — 1 and 



with L 3 (x) = - ln(l + x) + x — ^ 
has the following behavior: L^(x) 
L 3 (x) ~ -x 3 /3 for x < 1. 

Our main approximation will consist in replacement 
the function L 3 {x) by Loc(x), where — oo for x < — I 
and Lao = otherwise. The motivation for this approxi- 
mation is that we suppose that the main efffect of L 3 (x) is 
to forbid the region, where —ph(p) > 1. Then, minimiza- 
tion of the free energy functional amounts to satisfying 
the conditions 



ph{p) > -1 
h(r) > -1 
h(r) = — 1 for r < 1 



(3) 



which are in fact the minimum physical requirements for 
any correlation function h(r). In this sense we are build- 
ing a "minimum" model of a liquid. In addition to the 
constraints (|3|) we require that the function h{r) depends 
continuously on the density. The absence of a solution 



1 



which would be continuous in density would be a signal 
of a phase transition. This is, however, not found in the 
present calculations. 




1 2 3 4 5 



FIG. 1. Pair distribution function g(r) = h(r) + 1 of the 
3-dimensional model liquid, for densities (from top to bottom) 
p = 0.6, 0.7, 0.8, 0.9, 1.00, 1.1, 1.2, 1.23, 1.26, 1.29, 1.32, 1.35, 
1.38, 1.41, 1.44, 1.47, 1.5. The n-th curve from the bottom is 
shifted by n upwards. 

In 3D we can compute the function h(r) numerically by 
increasing slowly the density p and adjusting iteratively 
the function h(r) so as the conditions (|J) are satisfied. 
The resulting pair distribution function g(r) = h(r) + 1 
is shown in the Fig. The Fourier transform h[p) for 
p = 1.2 is shown in Fig. 0. 
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FIG. 2. Fourier transform of the correlation function for 
3-dimensional liquid, at the density p — 1.2. 

We can see that for densities up to about p = 1 the 
pair distribution function agrees qualitatively with the 
well-known results of HNC approximation or numerical 
simulations (sec U). However, at about p = 1.2 the 
behavior changes. A gap opens between the principal 
peak at r = 1 and the secondary peak at r ~ 2. the gap 
broadens with increased density and at about p ~ 1.5 



a second gap occurs around r ~ 2.2. We observed, that 
further compression leads to the occurrence of a third gap 
separating the peaks at r ~ 1.6 and r ~ 2. We expect 
that continued increase of the density results in increased 
number of gaps. 

The presence of the gaps is an artifact of the approx- 
imation. In reality the values of g(r) will not be strictly 
zero, but small. 

From the value of the radial distribution function at 
r = 1 the pressure can be computed Hand the result- 
ing equation of state is shown in Fig.T|, together with 
the results obtained by solving the HNC equations (|l|) 
and the formula computed in the Percus-Yevick (PY) 
approximation We can see that our model behaves 
qualitatively in the same way as the other approxima- 
tions, even though quantitative agreement is poor. On 
the other hand, the equation of state of our model does 
not differ from either HNC or PY approximation more 
than these two approximations differ one from the other. 
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FIG. 3. Equation of state for 3-dimensional liquid. Our 
model: full line (red). HNC approximation: dashed line 
(green). Percus-Yevick approximation: dotted line (black). 

We can see from these results, that the present ap- 
proach in 3D gives at least qualitatively sensible results. 
However, our aim is not to provide a new approximation 
for real three-dimensional liquids, but a model which de- 
scribes reasonably well the qualitative features of a liquid 
and is soluble in the limit of infinite spatial dimension. 
This will be done in the next section. 

III. SOLUTION OF THE MODEL IN HIGH 
DIMENSION 

In this section we will investigate a Z?-dimensional ver- 
sion of the model, with D ~ 2N + 3 and N — > oo. The 
main quantity of interest is again the correlation function 
h(r) = g(r) — 1. We will suppose that it has the form 
h(r) = ho(r) + h(r) where ho(r) = — 6(1 — r) and h(r) = 
for r < 1. 

The pressure is directly related to the value of that 
r = 1, more precisely to lim^^ h(r) = h(l + ) (sec Q) 



2 



(4) 



In the following, we will use rescaled quantities. We will 
use p = pVo for the density, while for the pressure we 
nse p = Y § P. 

In order to check the conditions (pft we should know the 
properties of the Fourier transform in high dimension. 
The relevant formulae used in this section are given in 
the Appendix A. 

In the zeroth approximation h(p) = —Vo ^&(p), which 
is correct as long as p < 1/Vb. So, for p < 1 we have 
h(r) — and the pressure is given by first virial correc- 
tion, 



For p > 1 we rewrite h in the form 

h(r)=h* (r) + h*(r) 



(5) 



(6) 



where still h* — h(r) for r > 1, but h*(r) is continuous for 
all r. The Fourier transform of Hq is easy to compute, if 
we know the value h*(l~) = — 1 — A, where A = h(l + ) = 
hi(l). Then 

h*(p) = -(l + A)V D * (p) ■ (7) 
In order to ensure h(p) > 1/p everywhere, we set 



h*(p) = 6i ( (A+l)V D * (p)-- 

P 



(8) 



(We denote 0\(x) = x for x > 0, 9\{x) — for x < 0.) 
The relevant quantity is 



Pi 



Inp 



(9) 



Indeed, one particle occupies space 2~ D Vn, so an abso- 
lute upper bound for the density is p\ < In 4. 

Let us suppose that p c < 1. This condition restrict 
the range of densities investigated to a certain interval, 
which will be found in what follows. 

For p < 1 the function 9 (p) — exp(— A^0 o (p)) is 
monotonously decreasing and we have the following equa- 
tion for p c . 



00 (Pc) = Pi ~ 



ln(^ + 1) 
N 



(10) 



When computing the inverse Fourier transform of 
h*(p) we need only the behavior around the point p c , 
which is 



h*(p)~-V D *' (p c )p c (l-?-) 

Pc 



(11) 



N\ VZp?+H-%(Pc))(A+l) 



D + l 

We obtain immediately the equation for A 



*o(p>) (12) 



2tt D + l 



V*p? +1 cb' (p c ))(A + l)^ (p c ). (13) 



In the interval p c < 1 we have 

Pc 



4>'q{pc 



< Pc 



(14) 



Hence, if we suppose that a solution such that A < 1 
exists, we have 



A < 



V£ft + W (p c ) * eM-2NK(p c )) (15) 



where 



K(p) = ln(l + ^Jl-p 2 ) - y/T 



hip 



(16) 



We have K(l) = and K(p) = -y/l -p 2 jp < 0, so 
K(p) > for p < 1. For fixed p c < 1 and N — > 00 
we have A < 1 and therefore we can neglect A in the 
equations ([To]) and (|l3|). 

So, we can conclude that in the range of densities p\ < 
Pic = ^o(l) = 1— ln(2) = 0.3068... the following equation 
for p c holds 



Pi 



= ln( 1 + s]\ -pi) - y/l-fi + 1 - In 2 



(17) 



The solution of the latter equation can be easily ob- 
tained in the form of power series. We show here only 
first several terms. The expansion up to order 16 is given 
in the Appendix B and the graph is shown in Fig. ^. 

f c = 4 Pl - 2 Pl 2 - \ Pl 3 - \ Pl 4 - 0( Pl 5 ) . (18) 
3 6 




FIG. 4. Dependence of the momentum p c on the rescaled 
density in the regime p c < 1. 



and gives, using (37) 



3 



From the solution Eq. (n~7]) we can compute the pres- 
sure. We write P — p + \jr + P c and re-scale the correc- 
tion as P c — exp(NPi). We obtain 



Pi = 2 - 21n2 + lnp2 



(19) 



We can see that the density p u = <f> (-) = 0.14676... 
separates two regimes. For p\ < p\ t the correction P c to 
the lowest virial becomes negligible for large TV, while for 
Pi > pit the correction diverges for TV — > oo. The graph 
of the function Pi (pi) is shown in Fig. 




FIG. 5. Equation of state for our model in the limit of 
infinite dimension, in the range of densities where p c < 1. 

For the correlation function in the interval rp c < 1 we 
have h(r) = exp(Nhi(r)), where 



hi(r) = 

l-21n2 + ln^-pi 



(20) 



+yl~— p 2 r 2 — ln(l + \/T — p 2 r 2 ) 



and for rp c > 1 we can use the following scaling h(r) = 
exp(Nhi(r)) cos(iV/i2(r)) where, using the expressions 
(B2n and (p3), we obtain 



(21) 



hi (r) = 1 — 2 In 2 + In p c — p\ — In r 



^{f) = y/Pc 1 " 2 ~ 1 — arctan yp%r^ — 1 



(22) 



We can also see that \h\(r) <c 1 for p\ < p lc , so that 
g(r) > for r > 1 and the gaps in <?(r), discussed in the 
last section do not occur. However, when p\ approaches 
pic the absolute value of h\{r) can be of order O(l) and 
a gap can appear at the density p\ c . The detailed inves- 
tigation of this process and the behavior of the model for 
Pi > pic remains an open question. 



IV. CONCLUSIONS 

We investigated a simple model for a hard-sphere liq- 
uid. By numerical solution in 3 dimensions, we found a 



qualitatively realistic behavior. The results for the equa- 
tion of state are compatible with the hypernetted chain 
and Percus-Yevick approximations. 

We solved the model analytically in the limit of large 
spatial dimension. We found that two scales of density 
and pressure appear, which corresponds to two regimes 
of density. For p < 1 the equation of state is given by first 
virial correction (||) , while for p > 1 the quantity relevant 
to the further virial corrections is p\ = In p/N and the 
pressure correction itself scale as P± = In P c /N (JTg|) . We 
have found the solution in the interval < Opi < p\ c = 
0.3068.... Two regimes are present within this interval. 
For pi < pu — 0.14676... the correction P c vanishes for 
large N, while for pi < p u it diverges for large N. 

It should be expected, that the presence of presence of 
gaps in g(r) will lead to qualitatively different behavior 
for densities higher than p\ c . 
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APPENDIX A 

Here we derive the formula for the Fourier transform 
in high dimension. 

Fourier transform in D = 2N + 3 dimensions is [x and 
k are D-dimensional vectors) 



f D (k)= / d D xe ixk f D (x) 



f D (x) = (2tt) 



—D 



d D ke~ 



(23) 



(24) 



We suppose that the functions depend only on the radial 
coordinate, f(r) = /d(x) for r = \x\ and f(p) = fo{k) 
for p = \k\. After rescaling p — p/N we will finally have 



M 



dr / dz 

i) 7-i 



where the constant Cn is fixed by condition that for 
f(r) = 9(1 — r) we have /(0) = Vn with Vd volume 
of the D-dimensional sphere, 



V D = 



2ir D / 2 
DT(§) 



e7r N N 
iV> 



(26) 



Similarly for the inverse Fourier transform we will have 
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m = 

/>oo /> 1 

Cn / dp dz 
with coefficient 



2(l+£) 



(1 



N 



N 



D 



m (27) 



(28) 



The calculation of the Fourier transform can be per- 
formed by the saddle-point method. The essential result 
is the Fourier transform of the surface of unit sphere, 
f(r) = S(r — 1). We obtain f(p) oc ^f{p) where 

*(P) = *o(p) = eM-NMP)) (29) 
for p < 1 and 

*(p) = = exp(-A^p))cos(iV<fe(p)) (30) 

for p > 1 . 

The explicit form of the functions <po,<f>i, <f>2 is given by 



6o(p) = 1 - In 2 + ln(l + ^1 - p 2 ) - ^1 - p 1 (31) 
6i(p) = 1 -ln2 + lnp (32) 

(33) 



4>2 (p) = V P 2 ~ 1 — arctan \ffP~- ~T 



Note that *(0) = 1. 

From here we can deduce the following Fourier trans- 
forms (8(x) — 1 for x > and = for x < 0). For 
/(r) = A 6(r - r): 



f(p) = AV D ftf V(pr ) 
For f(r) = A (l-r/r )6(r - r): 



(34) 



(35) 



Because the inverse Fourier transform has the same 
form and differs only in the factor C n instead of Cn , we 
can also immediately write for f(p) = A 9(p c — p): 



/M=(^) AV D p?*(p c r) 
and for f(p) = A (1 - p/p c ) 9(p c - p) 



AV D p L c 
D + l 



^>(p c r) 



(36) 



(37) 



vl = 

4 Pl 2 Pl 2 - 2/3 Pl 3 - 5/6 Pl 4 - | Pl 5 - ^| Pl 6 



6889 7 24721 8 



1260' 



2620169 9 64074901 10 
90720 91 ~ 907200 Pl 



1775623081 n _ 1571135527 12 
9979200 Pl 3421440 Pl 



1882140936521 13 70552399533589 14 
1556755200 Pl 21794572800 Pl 

_ 2874543652787689 15 25296960472510609 16 
326918592000 Pl 1046139494400 Pl ' 



(38) 



APPENDIX B 



[2. 

[3; 
[4; 

[5. 

[7; 

[8 
[9 
[10 

[11 
[12 

[13 
[14 

[is; 

[16 
[17 
[18 



Electronic address: slanina@fzu.cz 

C. A. Croxton, Liquid state physics - A statictical 
mechanical introduction (Cambridge University Press, 
1974). 

J. A. Barker and D. Henderson, Rev. Mod. Phys. 48, 587 
(1976). 

M. Mezard and G. Parisi, J. Phys. A: Math. Gen. 29, 
6515 (1996). 

M. Mezard and G. Parisi, Phys. Rev. Lett. 82 , 747 
(1999). 

M. Mezard and G. Parisi, J. Chem. Phys. Ill, 1076 
(1999). 

B. Coluzzi, M. Mezard, G. Parisi, and P. Verrocchio, J. 
Chem. Phys. Ill, 9039 (1999). 

B. Coluzzi, G. Parisi, and P. Verrocchio, J. Chem. Phys. 
112, 2933 (2000). 

M. Mezard, Physica A 2 56, 352 (1999). 
B. Doliwa and A. Heuer, |cond-mat/991242l| (1999). 
B. Coluzzi, G. Parisi, and P. Verocchio, Phys. Rev. Lett. 
84, 306 (2000). 



C. Dasgupta and O. T. Vails, |cond-mat/000208l| (2000). 
T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. Lett. 
58, 2091 (1987). 

T. R. Kirkpatrick and D. Thirumalai, Phys. Rev. B 36, 
5388 (1987). 

E. Marinari, G. Parisi and F. Ritort, J. Phys. A: Math. 
Gen. 27, 7615 (1994). 

E. Marinari, G. Parisi and F. Ritort, J. Phys. A: Math. 
Gen. 27, 7647 (1994). 

M. Mezard, G. Parisi, and M. A. Virasoro, Spin Glass 
Theory and Beyond (World Scientific, Singapore, 1987). 
A. P. Young (editor), Spin glasses and random fields 
(World Scientific, Singapore, 1997). 
E. Meeron, J. Math. Phys. 1, 192 (1960). 
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